Finite temperature excitations of Bose gases in anisotropic traps 
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The mode frequencies of a weakly interacting Bose gas in a magnetic trap are studied as a 
function of the anisotropy of the trap. As in earlier works the generalized Hartree-Fock-Bogoliubov 
equations within the Popov approximation (HFB-Popov) are used for our calculations. The new 
feature of our work is the combined use of a mode expansion in a finite basis and a semiclassical 
approximation of the highly excited states. The results are applied to check the accuracy of the 
recently suggested equivalent zero-temperature condensate (EZC) approximation which involves a 
much simpler model. 
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I. INTRODUCTION 



New experiments on condensed Bose gases with oscillating trap potential permit the excitation of low-lying 
collective modes with given symmetry. The measurement of their excitation frequencies provides a good opportu- 
nity to analyze the applicability of different approximations and numerical approaches in many-body-theory. Zero 
temperature calculations with different number of atoms in the condensate [|]|(| are in excellent agreement with the 
corresponding spectra measured in the JILA-TOP [0 and the MIT trap j|. The frequently used starting point for 
the theoretical studies are the coupled Bogoliubov equations. Their eigenvalues can be determined numerically by 
expanding the solutions in a basis-set of orthonormal functio ns ft[ $i ■ Analytical solutions of the Bogoliubov equations 
have been obtained in the hydrodynamic limit by Stringari ]10| and in ]TT] , p~2[ . 

For finite temperatures a possible and frequently used extension of the Bogoliubov equations are the Hartree-Fock- 
Bogoliubov equations with neglection of the anomalous expectation values, commonly referred to as Hartree-Fock- 
Bogoliubov (HFB)-Popov equations [p^|Jl4[ . Without their neglection the anomalous expectation values appearing in 
the self-consistent Hartree-Fock theory would lead to an unphysical energy gap of the collective modes in a spatially 
homogeneous Bose-condensed system violating the Hugenholtz-Pines theorem |{L5| . In a spatially inhomogeneous 
trapped Bose-condensate the Hugenholtz-Pines theorem does not apply, but in this case the Kohn-theorem jl6j would 
be violated, which states that the dipole excitations have the frequencies of the harmonic trap. 

The HFB-Popov equations have been solved within the local density approximation Jl^,^8|, which is sufficient to 
determine the thermodynamic properties of the trapped Bose-condensates. Beyond that, similar to the numerical 
method used in the case of the zero-temperature Bogoliubov equations, the excitation spectra can be derived by an 
expansion in a basis set of eigenfunctions. Dodd et al. [ [l9| present the numerical results together with the experimental 
data. For temperatures T < 0.65 T (T is the theoretical transition temperature for the ideal, trapped Bose gas) 
the calculated values differ by less than 5% from the experiment. On the other hand, it was also pointed out in Eg ] 
that the same accuracy between numerics and experiment is already obtained by a much simpler approximation. The 
excitation spectra at finite temperatures are simply derived by a zero temperature calculation but with the number of 
atoms in the trap reduced to account for the thermal depletion. Correspondingly, this solution is denoted as equivalent 
zero-temperature condensate (EZC)-solution [ |l9| |. 

In the present paper we present a further numerical study of the HFB-Popov equations with two aims: On the 
one hand we wish to calculate the excitation spectra by combining the use of a basis set expansion for the low-lying 
modes with the local density approximation for the high-lying modes beyond a suitably chosen energy cut-off. This 
eliminates truncation errors while keeping the discrete basis set of modes reasonably small, thereby speeding up our 
computation. The second aim then is to use this gain in efficiency for a more detailed check of the accuracy of the 
EZC-approximation in comparison with the HFB-Popov approximation. In particular we shall provide a comparison 
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for axially symmetric traps with arbitrary anisotropy both for low lying excitations and also for excitations at energies 
above the chemical potential. 



II. THEORY 



The HFB-Popov mean-field theory for inhomogeneous dilute Bose gases has been derived in detail by Griffin in Ref. 
|p^| . We merely give a brief summary of the basic equations. The Bose gases in the experiment are weakly interacting 
and dilute. Therefore, the s-wave approximation v(r — r') = g 5(r — r') is adequate to describe the interaction where 
g = A-KTi 2 a/m and a is the s-wave scattering length. To the Bose field operator ijj(r,t) the usual decomposition j2(J 
in a c-number part plus an operator with vanishing expectation value is applied: ip(r : t) = <&(r) + ip(r,t). $(r) 
represents the condensate wave function and the operator ifj(r,t) the excitations of the condensate. 

This ansatz is inserted in the Hciscnberg equation of motion (in units of Ti = 1): 



dip(r,t) ( V 
i — ' 



ill = (-^+ U ^p(r)-^(r,t)+g^(r,t)^(r,t)i>{r,t) (1) 

where {/trapC 7 *) = iri(UpP 2 + cu 2 z 2 )/2 is the trap potential (m the atomic mass and u> p and oj z the radial and axial trap 
frequencies) . The statistical average over (|l|) and the replacement of the cubic term ffi (r , t)ip{r, ^^(r, t) by the average 
in mean-field approximation 2 (ijj^(r, £)?/7(r,£)^) "0(r,t) neglecting the anomalous expectation value (?p(r, t)i/'(r, t)^ 
and its complex conjugate leads to the generalized Gross-Pitaevskii (GP) equation in the HFB-Popov approximation 

~ + CWr) - m) *(r) + 9 K(r) + 2n(r)] §(r) = . (2) 

Here we introduced the local density of the condensate n c (r) — |$(r)| 2 and of the depletion n(r) — (ip^(r, t) ip(r, t) 
The subtraction of equation (||) from (|l]) and a treatment of the cubic terms within the mean-field theory similar to 
the one of equation (|^) yields two coupled equations of motion for "0(r, t) and its adjoint. 

We only state the final results after the insertion of the Bogoliubov transformation V>(r,i) = 

J2j i u j( r ) &j e~~ lEjt + v*(r) ctj e lEjt ] leading to the generalized Bogoliubov equations for the generalized eigenfunctions 
Uj(r),Vj(r) and their corresponding eigenvalues Ej\ 

Cuj(r) +gn c {r) Vj {r) = E jUj {r) 

(3) 

t Vj (r) + g n c (r) Uj (r) = -E 3 v 3 (r) 
where £ = — + ?7trap(f) — /J, + 2 gn(r) with the total density n(r) = n c (r) + n(r) . 

The number density of particles fi(r) outside the condensate is given in terms of the thermal number of quasi-particles 



a) a-) = (e^/fer-i) 1 by 



Kr) = E{(M r )l 2 + l^( r )l 2 ) (4 & i) + l«iWI 2 } ( 3a ) 



where Uj(r), Vj(r) are normalized by J d 3 r ((^(r)! 2 — |wj(r)| 2 ) = 1. 

The coupled equations can either be solved by direct expansion of the solution in the basis set of the free trap potential 
as in |l9| or by using the local density approximation as in |l7],[l^]. The local density approximation is the leading 
order of a semiclassical approximation and amounts to setting 

Uj (r) -> u(p, r)e*M v 3 (r) - v(p, r)e^ ^ . . . / . . . . (4) 

where <j> is defined by V</> = p and u(p, r), v(p, r) are normalized by 
KP,r-)| 2 -Kp,r)| 2 = l. 

In the semiclassical limit the functions u(p, r) and v(p, r) are slowly varying on the scale of the harmonic oscillator 
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length d — {l/muj) 1 / 2 with lo — (w^z) 1 / 3 , and the derivatives of u and v as well as the second derivatives of <f> are 
negligible. The equations in (||) are then reduced to the algebraic form 



— + ^trap(r) — fi + 2gn(r) ) u(p, r) + gn c (r)v(p, r) = e(p, r)u(jp, r) 



P 2 

Y~ + ^trapO) - /i + 2gn(r) ) v(p, r) + gn c {r)u{p, r) = -e(p, r)v(p, r) . 



(5) 



The solution is straight-forward and we are led to the excitation spectrum e(p, r) — {e 2 HF {p, r) — g 2 rt^(r)) 1 / 2 with 

2 

the Hartree-Fock energy 6hf{p, t) = f— + U trap (r) — fj, + 2gn(r). The non-condensate density is then obtained by a 
simple integration over the momenta: 



d 3 p \e H F(p, r) ( 1 1\ 1 



(2tt) 3 



e(p,r) \exp{e(p,r)/k B T)-l 2 J 2 



e(e^ F (p,r)-^(r)) . (6) 



Let us first consider the case where the local density approximation is used for the whole excitation spectrum, not only 
the high- lying part. Then it is necessary for consistency to treat also the condensate in the corresponding approxi- 
mation, which is the finite-temperature Thomas-Fermi approximation. It applies to sufficiently large condensates. In 
this case the kinetic energy term in eq.(||) is neglected and one obtains with n(r) = h(r) + n c (r) 

Mr) - (M - t/tra P (r) - 2gh(r)) Q(fi- U trap (r) - 2gh{r)) (7) 

and hence e(p, r) reduces to 

{/ 2 X 1 / 2 

[fjk + 9 n c(r)) - g 2 n 2 c (r)J if M > U tvap (r) + 2gn(r) (g) 
^7 + U tTap (r) - /i + 2gn{r) if fi < U tTa , p (r) + 2gn{r) 

Inside the condensate the excitation energy is then gapless like in the spatially homogeneous case. Moreover the 
integrand in (Q) for r inside the condensate and p — > becomes kBTm/p 2 and therefore the integral (^) converges 
for small p (besides, of course, converging also for large p). 

On the other hand, even if the Thomas-Fermi approximation for the condensate is not applicable because the 
condensate is to small one may still apply the semiclassical local density approximation for sufficiently high lying 
states. But for the low-lying states the local density approximation is then inconsistent, as can e.g. be seen from the 

2 2 1/2 

fact that e(p, r) inside the condensate now has a gap for p — > of size E gap (r) = ^ ^ m ^j | j m ^^j + 2<?n c (r)^ 

in space points r where V 2< I , (r) > 0. Therefore, one has to introduce an energy cut-off e c in this case, below which 
the excitations are treated by solving eqs.(|J) exactly, expanding in a suitable basis, and above which the local density 
approximation can still be employed. Then eq.(^) is replaced by 

den(e,r) (9) 



where 
and 

n(e, r) 



Sj(r) = Mr£^ + K(r)l * (10) 



l 3 / 2 [ 1 1 



V2n 2 \ W fc BT _ i 2 2^e 2 + g 2 n 2 {r) j 
x \J \/e 2 + g 2 n 2 (r) - U tlap {r) + fi - 2gn(r) 



ill) 



It is our goal here to implement this hybrid procedure. The choice of the energy cut-off depends on the purpose 
of the calculation. If the aim is the calculation of thermodynamic properties requiring only the density of states it is 
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sufficient to take the energy cut-off low at an energy of a few w, and to restrict the discrete part of the sum over states 
to only very few discrete modes. If on the other hand the goal is to calculate the mode spectrum of the low-lying 
modes accurately up to a given energy, then this energy determines the cut-off. However, the employed basis set must 
then still be chosen considerably larger (in practice about twice as large) than the number of discrete modes to be 
calculated, in order to make sure that no discrete mode below the cut-off is missed. In this way it is ensured that 
eq. (||) really gives an accurate value of the particle density outside the condensate. The advantage of our hybrid 
approximation over the hitherto used simple cut-off of the employed basis set is that the truncation error in h(r) is 
avoided, which allows us the use of a smaller basis set and therefore a more economical computation. 



III. COMPUTATION 

The self-consistent procedure we use combines and extends previously used methods of ref. |g0J3- For a gi ven 
number of atoms, temperature and anisotropy of the trap we start e.g. with h = or a previously determined better 
estimate for n, solve eq.(||) for $(r), n c (r) — \$(r)\ 2 and fi from the condition that N — J d 3 r (n c (r) + u(r)), solve 
then eq.(|^) for all states up to the chosen energy cut-off to find Uj(r), Vj(r) and Ej, and determine finally an improved 
value for n(r) from eqs.(||), (|l0|). 

To calculate the second term in eq.(^|) it is necessary to introduce a finite region in space, where we perform the 
energy integral at certain points. Due to the axial symmetry of our problem this finite region is a rectangle in the r 
and z plane. This rectangle is chosen in such a way that n(r) and n c {r) fall fast enough to zero approaching the edges. 
In this rectangle we introduce a fine grid and the energy integrals are performed at every grid point numerically. Since 
both, n c (r) and n(r) are smoothly varying functions the grid methods are suitable from the numerical point of view. 
We use the same grid in solving the Gross-Pitaevskii eq.(||). To reach fast convergence for the solution of (||) we apply 
the Newton-Raphson algorithm in combination with multigrid methods as discussed in ref. [^2| . 

However, the grid methods are not appropriate in determining the solutions of eq.(|^) especially for the high- lying 
states. For this reason we solve the Bogoliubov-equations in the basis of the eigenfunctions of the free trap. In this 
basis we only need to calculate numerically the matrix elements of gn c {r), gn(r) in harmonic oscillator eigenstates. 
These matrix elements are decaying fast enough, because these quantities are smooth and localized around the origin. 
For our purposes the best method to solve (||) in our basis is that of Hutchinson et al. H . Then we do not have to 
introduce a basis for Uj and Vj separately, a single basis set is enough but two subsequent diagonalizations must be 
performed. For further details see Rj. 



IV. RESULTS 



We shall now discuss some results obtained by the procedure described in the previous section. In the calculations 
we present here we choose N = 2000 Rubidium atoms and T = 0.5 T c where T c is the critical temperature of an ideal 
Bose-gas with the same parameters. The anisotropy of the trap is described by the parameter 



i^j2 I o 

P = 2 7T^-, %■> U p = V av \l- — , U z = Wav V 1 + P (12) 

2uj z p + V 2 

where uj av is fixed to to av = 137 x 27r Hz. As a check on our code we reproduced and confirm the zero-temperature 
results for the low- lying excitation spectrum as a function of the anisotropy parameter presented in fig. 1 of ref. ||. 
A further check on our code, now at finite temperature, is the agreement of our results for the low-lying excitation 
spectrum as a function of temperature for the special anisotropy of the JILA-TOP trap with the results in [[l9| . 

Let us now turn to our own results. In fig. ^ we fix j3 at [3 — 1.4 and plot the spectral distribution g{e) of the 
particles in the trap over energy e in units of ui p . Thus AN = g(e)de is the number of atoms with energy e in the 
interval de. In fig. [j] we compare the spectral distribution 

9ld(e) = J d 3 rfi(e,r)e{e~U t r ap (r)+n-2gn(r)) (13a) 

obtained in the local density approximation with n(e,r) from eq.(|TI|) (full line), with the spectral distribution 

(13b) 



g d (e) ^J2f dVfi^r)— !=««? f-(^) 2 ) 
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obtained from the low-lying discrete modes by folding n, (r) ( |ic| ) with a Gaussian of variance U) p /\f2. Some smoothing 
on the quantum scale is necessary for a meaningfull comparison. The energy cut-off was chosen at e c = 30 lu Pi which 
is the reason for the fall-off of gd at around this energy. It can be seen that the local density approximation is rather 
good for energies somewhat above the chemical potential and higher. For energies at fi or lower the local density 
approximation increasingly fails and the spectral distribution from the discrete modes must be used. The oscillations 
by which <?d(e) differs from <?id(e) may be understood semiclassically f23[| by analysing the short periodic orbits of the 
classical dynamics generated by the Hamiltonian e(p,r) p4| . 

Next we discuss results for the excitation spectrum at temperature T — 0.5 T c as a function of the anisotropy 
parameter (3. We use the same parameter values as chosen in ref. [|l9| where the excitation spectrum for the special 
value of (3 = 1.4 of the JILA-TOP trap was presented. This case is therefore contained in our results. In fig. || we 
show the behavior of the chemical potential /i for N — 2000 atoms as a function of (3. In the vicinity of the limiting 
cases (3 = —1,2 the chemical potential changes rapidly and is expected to tend towards the groundstate energies of 
the free trap. 

The parity quantum number and the azimuthal quantum number m are good quantum numbers for axially sym- 
metric harmonic traps. On figs. ^-^| we show the elementary excitations as a function of (3 for different parity classes 
and m. For clarity only the excitation spectra with the quantum numbers m = 0, 1, 2, 3 are presented. On Figures ^ 
and ^ the lowest curves show the behavior of the center of mass modes, for odd parity with m — and for even parity 
with m = 1. They represent the collective oscillations of the atoms with the harmonic frequencies of the free trap 
ui z and Co p . By Kohn's theorem fl6f| they must occur in all finite temperature spectra and would be straight lines 
in our pictures. From our numerical data we see that Kohn's theorem is not exactly maintained in the HFB-Popov 
approximation, but the discrepancy for all anisotropies is so small, that it is practically negligible at T — 0.5 T c . In 
the excitation spectra (figs. |^-^|) modes with different axial quantum numbers m cross freely as one would expect if (|^) 
is considered as a linear set of eigenvalue equations. However, due to the coupling of these equations via n(r), which 
depends on the |itj(»")| 2 and the |uj(r)| 2 , the linearity is actually broken for finite temperature. But the nonlinearity 
enters only via a global coupling. Therefore, in discussing the crossings of levels we may consider n(r) and n c (r) 
as given by their final self-consistent values in which case (||) does become linear explaining the free level-crossings 
for different quantum numbers m. Some avoided crossings between levels with the same m can also be seen in these 
figures. 

A remarkable feature of the above spectra is that the levels tend to a few common eigenvalues in the limits (3^2 
or (3 — > — 1. In the first case (J3 — 2) uj p goes to zero and the limiting energy eigenvalues do not depend on the 
azimuthal quantum number m but still on the parity. In the opposite case ((3 — —1) they now depend on m but not 
on the parity. However, in the limit [3 — 2 the spectra agree with that of the free harmonic oscillator with u z = 0. 
This fact can be easily understood because both, the condensate and the thermal density vanish due to the repulsion 
between particles and the vanishing confinement in radial direction. The same argument can be applied for (3 = — 1 
with vanishing confinement in the z-direction. However, in the close neighborhood of (3 = — 1 some quantities (see 
for example the behavior of the chemical potential in fig. |J) change too rapidly and it is extremely difficult to get 
reliable data for the spectra in that region. Note that the limiting situations are quite different from those studied 
by Stringari p5fl (see also ref. Jl2|) under the conditions that the hydrodynamic approach is applicable, though some 
features are common. 

The low-lying excitation spectrum as a function of the anisotropy parameter for zero temperature has been presented 
by Hutchinson et al . The results are remarkably similar to ours considering the widely different physical conditions 
under which they are obtained. The similarity of the results of finite temperature calculations in the HFB-Popov 
approximation to zero-temperature results was first noted by Dodd et al |19|| . For anisotropy (3 = 1.4 they compared 
their results with an effective zero-temperature calculation (EZC), where N — N c , using as effective particle number 
N the number of particles in the actual condensate at temperature T . The close agreement they found led to the 
conclusion that in practice the excitation spectra are determined mainly by the atoms in the condensate, so that the 
thermal density h(r) can be neglec ted in the GP and the Bogoliubov equations (Q) and (|J). Furthermore, at least for 
the high- lying excitations |],|2(|g7]] Vj(r) can be set to zero in (H) @,^] and the equations (|J) reduce to an ordinary 
Schrodinger equation. We can now use our result for arbitrary values of (3 and looking also at higher lying excitation 
energies to provide a much more sensitive test of the accuracy of the EZC model. 

In fig. H we present the low-lying EZC mode spectra for the same case as in fig. ||. In the EZC approximation 
the Kohn-theorem is satisfied exactly, as there is no thermal cloud in the effective zero-temperature system. The 
excellent agreement with the HFB-Popov results, even in fine details, over the whole interval of anisotropies confirms 
the remarkable usefulness of the EZC approximation for the calculation of the low-lying excitation spectrum. 

However, considering a broader range of energies one can find differences between the predictions of the HFB-Popov 
and the EZC-calculations. In order to make a quantitative comparison we fix the anisotropy parameter (3 in fig. u\ to 
the value of the JILA trap j3 — 1.4. Here we plot the energy-differences between the HFBP- and the EZC-calculations 
Aei = chfbp ~ £ezc and between the HFBP- and the Hartree-Fock-calculations Ae2 = chfbp — £hf as a function of 
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the HFBP-excitation energies. At very low energies Aei is very small confirming again that in the low-lying spectra 
EZC levels agree with the HFB-Popov levels. With increasing energy the EZC reveals an increasing shift and Aei 
shows larger fluctuations. The difference between the eigenvalues is less than 0.2oj av . It is not too big, but shows 
clearly the magnitude of the difference between the two calculations. For energies above twice the chemical potential 
the Hartree-Fock model (vi neglected, but including the thermal density n) turns out to be a better approximation 
for the HFB-Popov equations. The shift for large energies can be understood in the following way: There is a shift in 
the chemical potential A[i = /^hfb — Mezc = 0.19 because of the neglection of the almost constant thermal density in 
the condensate region in the EZC. Considering only small energies the eigenmodes are localized inside the condensate 
where the term 2gn and the chemical potential shift cancel each other. Therefore the energy difference Aei is small 
for low energies. But for states localized outside, where even h is negligible the difference A^i is the main reason for 
the shift of the levels in the EZC calculations. The Fig. |^ is in accord with the expectation that the assymptotic limit 
of the energy shift is A/i. 



V. CONCLUSIONS 



We have investigated the Hartree-Fock-Bogoliubov-Popov equations to obtain the finite temperature excitation 
spectrum of trapped condensed Bose gases. Our computational method combined the use of the local density approx- 
imation above a certain cut-off and a solution of the discrete Bogoliubov equations below the energy cut-off. We have 
explored in a systematic way the excitation frequencies within a large range of parameters for the anisotropy of the 
trap potential. We gave numerical evidence for degeneracies occuring at extreme anisotropic situations (3 = 2, — 1 and 
argued that these two spectra can be explained by the non-interacting model. Further we compared the spectra of 
the HFB-Popov with the EZC calculations and showed that for large energies the EZC spectra are shifted by a small 
amount which is mainly due to the shift in the chemical potentials. 
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FIG . 1. The spectral distribution gid(e) ([l3a|) in the local density approximation is compared with the spectral distribution 



<7d(e) (13b) obtained by the expansion in a basis set for (3 = 1.4, TV = 2000, T = 0.5 T c . 

FIG. 2. fi in units of ui a v as a function of the anisotropy parameter for N = 2000, T = 0.5T C . 

FIG. 3. The lowest HFB-Popov excitations lo 2 in units of uj„ as a function of (3 for the even-parity modes with m = (open 
circles), m = 2 (black stars). Other parameters are N = 2000, T = 0.5 T c . 

FIG. 4. The same as in fig. |^ for the even-parity modes m = 1, m = 3. 

FIG. 5. The same as in fig. ^ for the odd-parity modes m = 0, m = 2. 

FIG. 6. The lowest equivalent zero condensate excitations u 2 corresponding to fig. ^| in units of ui 2 v as a function of j3 for 
the even-parity modes with m = (open circles), m — 2 (black stars). 

FIG. 7. The energy differences Aei = ehfbp — £ezc (open circles) and Ae2 = ehfbp — £hp (crosses) as a function of the 
HFBP-energies. 
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